Snap ATAC data processing

  library(ggpubr)

Loading required package: magrittr

Attaching package: ‘magrittr’

The following object is masked from ‘package:purrr’:

    set_names

The following object is masked from ‘package:tidyr’:

    extract


Attaching package: ‘ggpubr’

The following object is masked from ‘package:cowplot’:

    get_legend

f74.rna.seu@meta.data <- f74.rna.seu@meta.data %>%
  rownames_to_column("cell") %>%
  left_join(annotation.df) %>%
  column_to_rownames("cell")

Joining, by = "cell"

First I have made the snap file from the cellranger fragments.tsv file. Script is multiOmic_benchmark/preprocess/fragments2snap.sh.

Following integration vignette

Barcode selection

Filtering based on number of reads per cell and ratio of fragments that are within promoters

barcode.ls = lapply(seq(snap.files), function(i){
  barcodes = barcode.ls[[i]];
  idx = which(
      barcodes$logUMI >= cutoff.logUMI.low[i] & 
      barcodes$logUMI <= cutoff.logUMI.high[i] & 
      barcodes$promoter_ratio >= cutoff.FRIP.low[i] &
      barcodes$promoter_ratio <= cutoff.FRIP.high[i]
  );
  barcodes[idx,]
});
x.sp.ls = lapply(seq(snap.files), function(i){
  barcodes = barcode.ls[[i]];
  x.sp = x.sp.ls[[i]];
  barcode.shared = intersect(x.sp@barcode, barcodes$barcode);
  x.sp = x.sp[match(barcode.shared, x.sp@barcode),];
  barcodes = barcodes[match(barcode.shared, barcodes$barcode),];
  x.sp@metaData = barcodes;
  x.sp
})
names(x.sp.ls) = sample.names;
x.sp = Reduce(snapRbind, x.sp.ls);
x.sp@metaData["sample"] = x.sp@sample;
x.sp

number of barcodes: 6043
number of bins: 0
number of genes: 0
number of peaks: 0
number of motifs: 0

Add cell-by-bin matrix

x.sp = addBmatToSnap(x.sp, bin.size = 5000)

Epoch: reading cell-bin count matrix session ...

Binarize matrix

Some items in the count matrix have abnormally high coverage perhaps due to the alignment errors. Therefore, we next remove top 0.1% items in the count matrix and then convert the remaining non-zero values to 1.

x.sp = makeBinary(x.sp, mat="bmat")
x.sp

number of barcodes: 6043
number of bins: 651791
number of genes: 0
number of peaks: 0
number of motifs: 0

Filter bins

Filter out bins overlapping w ENCODE blacklist

black_list = read.table("~/annotations/hg38.blacklist.bed.gz")
black_list.gr = GRanges(
  black_list[,1], 
  IRanges(black_list[,2], black_list[,3])
);
idy = queryHits(
  findOverlaps(x.sp@feature, black_list.gr)
);
if(length(idy) > 0){
  x.sp = x.sp[,-idy, mat="bmat"];
};
x.sp

number of barcodes: 6043
number of bins: 651765
number of genes: 0
number of peaks: 0
number of motifs: 0

Exclude bad chromosomes

chr.exclude = seqlevels(x.sp@feature)[grep("random|chrM", seqlevels(x.sp@feature))]
idy = grep(paste(chr.exclude, collapse="|"), x.sp@feature)
if(length(idy) > 0){
  x.sp = x.sp[,-idy, mat="bmat"]
}
x.sp

number of barcodes: 6043
number of bins: 650342
number of genes: 0
number of peaks: 0
number of motifs: 0

remove the top 5% bins that overlap with invariant features such as the house keeping gene promoters

bin.cov = log10(Matrix::colSums(x.sp@bmat)+1)
# bin.cov = Matrix::colSums(x.sp@bmat)
hist(
  bin.cov[bin.cov > 0], 
  xlab="log10(bin cov)", 
  main="log10(Bin Cov)", 
  col="lightblue", 
  # xlim=c(0, 5),
  breaks=100
);

bin.cutoff = quantile(bin.cov[bin.cov > 0], 0.95)
idy = which(bin.cov <= bin.cutoff & bin.cov > 0)
x.sp = x.sp[, idy, mat="bmat"];
x.sp

number of barcodes: 6043
number of bins: 539584
number of genes: 0
number of peaks: 0
number of motifs: 0

Remove any cells of bin coverage less than 1,000. The rational behind this is that some cells may have high number of unique fragments but end up with low bin coverage after filtering. This step is optional but highly recommanded.

idx = which(Matrix::rowSums(x.sp@bmat) > 1000);
x.sp = x.sp[idx,];
x.sp

number of barcodes: 5793
number of bins: 539584
number of genes: 0
number of peaks: 0
number of motifs: 0

Dimensionality reduction

Uses diffusion map algorithm w sampling technique to make it fast.

## Sample 100 cells as landmarks 
row.covs.dens <- density(
  x = x.sp@metaData[,"logUMI"], 
  bw = 'nrd', adjust = 1
)
sampling_prob <- 1 / (approx(x = row.covs.dens$x, y = row.covs.dens$y, xout = x.sp@metaData[,"logUMI"])$y + .Machine$double.eps);
set.seed(1)
idx.landmark.ds <- sort(sample(x = seq(nrow(x.sp)), size = 1000, prob = sampling_prob))
## Split between landmark and query cells
x.landmark.sp = x.sp[idx.landmark.ds,];
x.query.sp = x.sp[-idx.landmark.ds,];
## Run diffusion map on landmark
x.landmark.sp = runDiffusionMaps(
  obj= x.landmark.sp,
  input.mat="bmat", 
  num.eigs=50
);

Epoch: checking the inputs ...
Epoch: computing jaccard similarity matrix ...
Epoch: fitting regression model ...
Epoch: performing normalization ...
Epoch: computing eigen decomposition ...
Epoch: Done

x.landmark.sp@metaData$landmark = 1;
## Project query cells
x.query.sp = runDiffusionMapsExtension(
  obj1=x.landmark.sp, 
  obj2=x.query.sp,
  input.mat="bmat"
)

Epoch: checking the inputs ...
Epoch: computing jaccard similarity matrix ...
Epoch: performing normalization ...
Epoch: projecting query cells to the reference ...
Epoch: Done

x.query.sp@metaData$landmark = 0;
## Combine 
x.sp = snapRbind(x.landmark.sp, x.query.sp);

'rBind' is deprecated.
 Since R version 3.2.0, base's rbind() should work fine with S4 objects

x.sp = x.sp[order(x.sp@metaData[,"sample"])]; #IMPORTANT

To determine significant diffusion components: > We use an ad hoc method by simply looking at a pairwise plot and select the number of eigen vectors that the scatter plot starts looking like a blob. In the below example, we choose the first 15 eigen vectors.

plotDimReductPW(
    obj=x.sp, 
    eigs.dims=1:50,
    point.size=0.3,
    point.color="grey",
    point.shape=19,
    point.alpha=0.6,
    down.sample=5000,
    pdf.file.name=NULL, 
    pdf.height=7, 
    pdf.width=7
  );

Clustering and visualization

x.sp=runCluster(
    obj=x.sp,
    tmp.folder=tempdir(),
    louvain.lib="leiden",
    seed.use=10,
    resolution=0.7
  );

Epoch: checking input parameters
Epoch: finding clusters using leiden

Visualization

Make gene matrix

saveRDS("~/my_data/cellranger-atac110_count_30439_WSSS8038360_GRCh38-1_1_0.snapATAC.RDS")

Error in saveRDS("~/my_data/cellranger-atac110_count_30439_WSSS8038360_GRCh38-1_1_0.snapATAC.RDS") : 
  'file' must be non-empty string

Repeat visualization on gmat

FeaturePlot(f74.seu, features = c('PTPRC','CD4','CD8A','CD8B','CD79A','FOXN1','EPCAM','PDGFRA','GNG4', 'FOXP3','RAG1','RAG2','NKG7','CCR7'), reduction = "umap.snap")

All cells have the same value (0) of CD79A.All cells have the same value (0) of RAG2.All cells have the same value (0) of NKG7.All cells have the same value (0) of CCR7.

DimHeatmap(f74.seu, dims = 1:6, balanced = TRUE,ncol = 1)

31 Oct 2019 09:28:40 [rsession-] CLIENT EXCEPTION (rsession-jovyan): (TypeError) : Cannot read property 'a' of null;
org/rstudio/studio/client/workbench/views/source/editors/text/ChunkPlotPage.java#52::execute
org/rstudio/studio/client/workbench/views/source/editors/text/ChunkPlotWidget.java#49::onBrowserEvent
com/google/gwt/user/client/DOM.java#1414::dispatchEvent
com/google/gwt/user/client/impl/DOMImplStandard.java#312::dispatchEvent
com/google/gwt/user/client/impl/DOMImplStandard.java#334::dispatchUnhandledEvent
com/google/gwt/core/client/impl/Impl.java#244::apply
com/google/gwt/core/client/impl/Impl.java#283::entry0
rstudio-0.js#-1::eval
Client-ID: 7beb3877-891d-4e16-98d5-f8cca1f72d5d
User-Agent: Mozilla/5.0 (Macintosh  Intel Mac OS X 10_14_6) AppleWebKit/537.36 (KHTML, like Gecko) Chrome/77.0.3865.90 Safari/537.36

Differential analysis

Expression of markers in RNA data

map(1:10, ~ plot_cluster_marker_ATACvsRNA(.x))

All cells have the same value (0) of MIR4300HG.All cells have the same value (0) of NOS1.All cells have the same value (0) of RPH3A.Cannot convert object of class listggarrange into a grob.All cells have the same value (0) of C5orf66-AS2.All cells have the same value (0) of CYP3A43.

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

Embedding on most variable genes in ATAC (not in gene expression)

Thoughts

• There seems to be a correspondance for gene expression and accessibility in many “accessibility markers” (see Cluster 6) –> maybe integration would work better if I don’t take only the most variable features in the RNA data.